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Abstract. In this work we present our first results of binary black hole circular 
orbits using COCAL, the Compact Object CALculator. Using the 3+1 decomposition 
five equations are being solved under the assumptions of conformal flatness and max- 
imal slicing. Excision is used and the appropriate apparent horizon boundary con- 
ditions are applied. The orbital velocity is determined by imposing a Schwarzschild 

■ behaviour at infinity. A sequence of equal mass black holes is obtained and its main 
_ physical characteristics are calculated. 
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1. Introduction 

I ■ One of the most important tests of Einstein's theory of general rela- 

tivity is the search for gravitational waves. A great effort both in the 
experimental and theoretical problems has been made and detection 
can happen almost any time. A highly probable scenario will be that the 
J> \ gravitational wave is coming from a binary system of two black holes 

or two neutron stars or a black hole/neutron star system. Therefore 
the extraction of a waveform that represents such configurations is an 
l/^ I important step towards detection. 

■ From the mathematical point of view assuming spacetime is foli- 
" ated by three dimensional hypersurfaces S^, Einstein's equations can 

be written as an initial value problem for the first and the second 
fundamental form of Sj. Then we get two sets of equations; one set 
that provides initial data and another that evolves them to acquire the 



m 

00 



^ ■ full spacetime. In ( Uryu and Tsokaros, 2012 ) we provided a method 



^ . to solve the former set of equations and here we elaborate on these 

solutions and identify those that represent circular orbits. Also we 
present some preliminary results regarding the physical characteristics 
of these solutions. 

The spacetime metric on is written in 3-|-l form as 

ds^ = g^^ydx^dx" = -a^dt^ + -iij{dx' + ^'dt){dx^ + pUt). (1) 
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We assume the spatial three metric 7ij on the sUce to be conformally 
flat = tp'^fij. Then the system to be solved, which are Hamilto- 
nian and momentum constraints and the spatial trace of the Einstein's 
equation, becomes 

= (2a) 

Aft = -2 a A^dj In ^ - ^didj^^ , (2b) 

A(aV) = Ja^'4-^''> (2c) 
o 

where A := did^ is a flat Laplacian. The fleld variables il^,a, and /3* 
are the conformal factor, lapse, and shift vector, respectively. We also 
assume maximal slicing to Sf , so that the trace K of the extrinsic cur- 
vature Kij = Aij + \^ijK vanishes. The conformally rescaled quantity 
Aij becomes 

4 = ^(^4 +5jft-^/^A/3'=) , (3) 

where the derivative di is associated with the flat metric /jj, and con- 
formally rescaled quantities with tilde are defined by Ai^ = Ai^ and 
/3* = whose indexes are lowered (raised) by fij (P^)- 



2. Overview of the algorithm 



In our previous paper (Uryu and Tsokaros, 2012), we presented the 



new code for computing equilibriums of astrophysical compact objects 
- COCAL, Compact Object CALculator. In the COCAL code we cover 
the initial hypersurface with spherical grids like the one that appears 
in FiglU Characteristic features are the inner spherical surface Sa of 
the patch with radius (small circle at the center of Fig{T]) , the outer 
spherical surface Si, with radius r;,, and and excised sphere Se with 
radius re- The role of Sa is to exclude the region near the black hole 
singularity. Boundary conditions must be provided there. The role of Sf, 
is to reach the asymptotic region. Boundary conditions are also imposed 
at Sf). Finaly Se is introduced to improve the angular resolution and 
reduce the number of multipoles for resolving the companion object. 
The boundary value at is copied from the sphere of the same radius 
as indicated in FiglU so that the equal mass binary black holes can be 
calculated. 

The method that we use to solve the partial differential equations is 



the Komatsu-Eriguchi-Hachisu (KEH) method (Komatsu, Eriguchi, and Hachisu, 1989) 
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Figure 1. The core COCAL coordinate system where the initial value equations are 
being solved. One dimension is supressed. The sphere at the center, surface Sa, 
corresponds to one compact object. The radius of coordinate patch doesn't reflect 
the actual size. 



which essentiahy uses the representation theorem with a suitable chosen 



kernel iteratively until a fixed point is obtained. In (Tsokaros and Uryu, 2007) 
the KEH method was adapted to handle multiple coordinate patches 
with appropriate boundary conditions. In COCAL the construction of 
the kernel is intimately related to the geometry of Fig{T]and the bound- 
ary conditions it satisfies. Denoting by Ba the ball or radius r^, Bi, the 
ball of radius r?,, and Bf, the ball of radius rg, it is Sa = dBa, Sb = dBb, 
Se = dBe. Our computational domain is V = Bf, — {Ba U i?e) and we 
have dV = SaU SbU Se- 

A typical boundary value problem (BVP) that we encounter is 

V2$ = S($) inV, C^ = f ondV (4) 

where $ can be any of the metric potentials and C a first order linear 



operator. Following (Jackson, 1975, chap 3) we write the solution as 



= x{x) + ^iNT(a;), (5) 

where 

$int(:e) = -^ I G{x,x')S{x')(fx' 

[ \G{x,x')V^^{x')-<^{x')VG{x,x')]dS'a. (6) 
47r Jqv 
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and xi^) is the corresponding homegeneous solution of the BVP 



V^X = in y, 



Cx = f — ^^INT on dV 



(7) 



In Eq. ([6]) G{x, x') is the flat Green's function that satisfles V^G(x, x') = 
— 47rJ(x — x'). Expanding G in multipoles on a spherical coordinate 
system we have 



X — X' 



£=0 m=0 

xP;"(cos0)P7"(cose')cosm((^ - if') 



(8) 



where the radial Green's function gi{r,r') is defined by 



9l{r, r') 



(9) 



with r> := max{r, r'}, r< := min{r, r'}, and the coefficients tm are 
equal to eo = 1, and = 2 for m > 1. 

2.1. Implementation for Robin-Dirichlet boundary 
conditions 

When one computes inversion-symmetric initial data or when enforces 
the inner surface Sa to be an apparent horizon a Robin type boundary 
condition for the conformal factor is obtained. The BVP that has to 
be solved is 



= in V 



dr 2r 



f 



where /, ^f, are known functions. The corresponding BVP for the ho- 
mogeneous solution is 



dx ^ X_ 
dr 2r 



[x] 



V^X = inV 
= f- 



9$INT _|_ ^INT 



r=rt 



dr 



2r 



(10a) 
(10b) 
(10c) 
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Since 7-^, r ^ ^ are the solutions of the radial part of the Laplacian, we 
write 



00 £ 



\\A(,mr ^ ^ + CimAcos{m(t)) + [Bimr ^ ^ + L>£m/] sin(m</>)] } (11) 

where Aim-, B^m, Cim, and Dim are constants. From boundary condi- 
tions Eq. (fTObl) . (fTncl) and using the orthogonality relations 



2lT 



we get 



sin^nKp) cos (m' (f))d(f) = 


cos{m(t)) cos{m (j3)d4> = — ^mm' 





Atmr^, ^ ^ + Qmr^ = (2£ + 1) x 

/ " - «>int),=„ PricosO) cos{m<^)d^ 
JO 



„ , „ I 2(2^ + 1) 

($5 - $int),=,, iT(cos e) sm{m4>)dn 



JO 



and 



-Aimr-^~'^ + Cemry = 2 



X 



TT /'27r 

Jo 
'0 Jo 



(9$INT ^*INT 1 r>mr n\ r j,\Jn 

/ — ) P^ [cosO) cos[m(p)dil 



-BimTa ^ ^ + Demri ^ = — X 

9$INT <1'INt\ /)\ • f A\jn 

J Pf (cos 8)sm(m(p) ail 

9r 2r J ^^^^ 
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When we solve the above system of equations with respect to A^^, 
Bimi Cirrii D^j^ and substitute back to Eq. ([TT]) we get 



oo £ 



X{x) 



— YY 

47r ^ ^ 



e^i^-^PricosS) X 



t=0 m=0 



+ m)! 

[ - $iNT)iT(cos e') cos[m{<J) - <j)')]dn' 



dr 2r J 



PPicos e') cos[m((/. - (t)')]dVL' \ .(12) 



where 



a^r) = (2£ + 



An) :.fer' 



bi{r) 



i+1 



rb 



f+1 / X 2e+l 



The final solution will be obtained from the iteration of 

= x{x) + ^iNT(a:) 

where $int is given by Eq. 



3. Binary black hole circular orbits 



To solve for BBH, Eq. ()2ap .()2b p . (j2c|) are supplemented with boundary 
conditions at infinity as 

^^1^ = 1.0, p\,^^ = 0.0, aUoo = l-0 (13) 
so as flat space time is acquired, and at the black hole excision surface 



Sa (Cook and Pfeiffer, 2004 Grandclement, 2010) with 



dr 2r 

f3' 



no. 



(14a) 

(14b) 
(14c) 
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Instead of the infinity, we impose the boundary condition (jl3p at the 
sphere Sb with the radius r = in the asymptotic region. The bound- 
ary condition Eq. ()14cp encodes the freedom to choose the initial shce. 
For BBH in quasiequihbrium the choice for the lapse is largely irrele- 
vant. We are taking no = 0.1 on Sa- In Eq. ()14ap is the unit normal to 
the sphere Sa- This equation enforces the sphere Sa to be an apparent 
horizon. Finally Eq. ()14bp ensures that the spheres are in equilibrium 
and also informs about the state of rotation of the BH. M* is the con- 
stant translational vector defined by (0, d, 0), where d is the coordinate 
distance between the center of mass of the binary system and the center 
of the black hole excision surface. (f)\ is the rotational vector with respect 
to the black hole. Parameter $7 corresponds to the orbital velocity of 
the system, and parameter Vtg to the local rotation rate of the BH 
that relates to the spin of the BH. For any values of 17, Vtg a solution 
is obtained and we have to choose those that correspond to circular 
orbits. In this preliminary work, parameter $7^ is set to zero, which 
approximately corresponds to irrotational (non spinning) BBH solu- 



tion. As it is discussed in ( Caudill, Cook, Grigsby, and Pfeiffer, 2004 ) 



the value of Q.s that leads to non spinning binaries can be found by 
solving Eq. (|2ap .()2bp.(|2c p iteratively until the quasi- local spin 



5^ = ^ / V''^ifc<Ai(,)d5', (15) 

vanishes. We will present the result with such adjustment of the spin 
in our forthcoming paper. 



Following ( Gourgoulhon, Grandclement, and Bonazzola I, 2002) the 



value of the orbital velocity for a circular orbit is obtained by requir- 
ing equality of the ADM mass and the Komar mass. In the conformally 
flat spacetime these are calculated from 

Madm = I di^dS\ (16) 

^71" Joo 

MKomar = ^ [ diadSK (17) 

In our code the surface integrals are performed over the sphere Sf,- 
Alternatively we can convert the surface integrals at infinity to volume 
integrals and surface integrals on the BH thus providing a consistency 
check to the accuracy of our solution. Using Eq. ([2a|) we get for the 
ADM mass also 

Madm = t^ / ^''AijA'^dV-^ [ dii^dS'-^ [ di^dS\ (18) 
167r Jv 27r J 27r J 
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The ADM mass given by Eq. (jl8p will be used to test the accuracy of 
solving the Hamiltonian constraint. 

The total angular momentum in a hypersurface St is defined as 

J=^ f (K^- K6iM^dS, = ^f A]^^dS,. (19) 



oo 

,3 



Again this integral in our code is taken over the sphere Sb- 4>cm is 
the rotational vector with respect to the center of mass. Similarly we 
convert the integral at infinity to an integral at the throat plus a volume 
integral. The latter vanishes due to the momentum constraint Eq. ()2bp 
and the fact that (/)cm is a Killing vector for the flat metric. Finally we 
have 

J = -^ [ ^''^)^cmdS^ [ ^l^^'A^cl^l^dS,. (20) 

Comparing Eq. (j20|) with Eq. (jl9p we get an estimate for the violation of 
the momentum constraint. As a final check of the overall computation 
we consider the Smarr formula 

Madm - 217 J = I il^^diodS' [ ^^diadS\ (21) 

that relates the ADM mass, the angular momentum, and the computed 
orbital velocity. Typically by inserting Madm from Eq. (|16p and the 
calculated orbital velocity we can get a third value for the angular 
momentum of the system that can be compared against Eq. ()20p . and 
Eq. (USD. 

Another important quantity is the irreducible mass and the binding 
energy of the system. It is M\^^ = mi + m2 with rrii = \J Ai/Wir and 

Ai= [ ^^dS. (22) 



^irr- 



The binding energy is then = Madm — Afii 

Using the H3 resolution in Table [H a sequence of equal mass black 
holes is obtained whose main characteristics (separation parameter, an- 
gular velocity, ADM mass, binding energy, and angular momentum) are 
shown in Table |TI1 The renormalization was done using the irreducible 
mass. For the BH coordinate separation we used dg = 2.5. 

For two point particles of individual mass m moving in circular 
orbit of radius ii, Kepler's third law gives 0^ = 2m/{2R)^ , and since 
the total angular momentum of the system is J = 2mQR'^, we have 
that in Newtonian mechanics the total mass of the system M = 2m, 
the total angular momentum J, and the angular velocity J7 satisfy 
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Figure 2. The quantity 4J(r2/M^)^''^ is plotted against the separation parameter 
D. 



The deviation from the Newtonian value of the quantity on the left 
side of Eq. (j23p can be seen in Fig. [21 The closer the black holes are, 
the larger the difference between the Newtonian and general relativistic 
prediction. 

Plots of the ADM mass and the angular momentum versus the angu- 
lar velocity can be seen in Fig. El As the black holes come together the 
mass and the angular momentum exhibits a minimum that signifies 
the innermost stable circular orbit (ISCO). As we observe from the 
top panel of Fig. [3l the minimum of the ADM mass is different from 
the minimum of the angular momentum (middle panel). The reason 
for this discrepancy is that the magnitude of the spin as computed by 
Eq. (jlSp is not exactly zero and therefore the sequence is not strictly 
speaking irrotational. This is the reason that the characteristic cusp in 
the mass versus angular momentum plot is absent (bottom panel). As 
we discussed above this issue is resolved when we iterate over Vtg so as 
to make the spin, Eq. (fT5]) . to be zero. 



4. Discussion 

We have successfully computed a sequence of conformally flat initial 
data for non-spinning equal mass BBH solutions in circular orbits. 
Several authors have calculated sequences of this kind as models of 
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Figure 3. Top panel: plot of ADM mass versus the orbital velocity. Middle panel: 
plot of angular momentum versus the orbital velocity. Bottom panel: plot of ADM 
mass versus angular momemtum. 
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Table I. Grid parameters used in the computation of equal mass BBH 
sequence. {Nr, Ng, N,i,) are the numbers of grid points in {r,0,<j)) co- 
ordinates, and L is the number of multipoles included in the Green 
function. We vary the radius Ta in the range shown in the table to 
compute the solution sequence from larger to smaller separations. 



Type ra 






ds 


Nr 


Ne 




L 


H3 0.1-0.23 




1.125 1.25 


2.5 


384 


96 


96 


12 



Table II. Sequence of equal mass black holes on a maximal slice. The 
normalization is done using the irreducible mass. 



D 




MADM/Mi„ 




J/M2„ 


25.00 


0.04516 


0.98367 


-0.01633 


0.84973 


22.73 


0.05175 


0.98274 


-0.01726 


0.82913 


20.83 


0.05856 


0.98190 


-0.01810 


0.81218 


19.23 


0.06558 


0.98116 


-0.01884 


0.79820 


17.86 


0.07281 


0.98054 


-0.01946 


0.78667 


16.67 


0.08023 


0.98004 


-0.01996 


0.77721 


15.63 


0.08783 


0.97968 


-0.02032 


0.76953 


14.71 


0.09563 


0.97948 


-0.02052 


0.76343 


13.89 


0.10360 


0.97945 


-0.02055 


0.75874 


13.16 


0.11178 


0.97963 


-0.02037 


0.75538 


12.50 


0.12015 


0.98005 


-0.01995 


0.75327 


11.90 


0.12874 


0.98074 


-0.01926 


0.75240 


11.36 


0.13757 


0.98176 


-0.01824 


0.75281 


10.87 


0.14669 


0.98320 


-0.01680 


0.75460 



BBH inspiral due to the graviational wave radiation (GGB II, 2002 



Caudill, Cook, Grigsby, and Pfeiffer, 2004: Grandclement, 2010). They 



used spectral methods in their computations, and produced numerical 
solutions in higher precision compared to ours. In the COCAL code, we 
use standard, mostly second order, finite difference scheme. Our method 
is much simpler than the spectral method, and hence it is easier to ex- 
tend our code to include magnetic fields or neutron stars. Also we have 
demonstrated that the solutions are accurate enough to reproduce the 



results of (GGB II, 2002 Caudill, Cook, Grigsby, and Pfeiffer, 2004 Grandclement, 2010) 



Further details of the COCAL code will be discussed elsewhere. 
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